Density functional theory for fermions close to the unitary regime 
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We consider interacting Fermi systems close to the unitary regime and compute the corrections 
to the energy density that are due to a large scattering length and a small effective range. Our 
approach exploits the universality of the density functional and determines the corrections from the 
analyical results for the harmonically trapped two-body system. The corrections due to the finite 
scattering length compare well with the result of Monte Carlo simulations. We also apply our results 
to symmetric neutron matter. 
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Ultracold fermionic atom gases have attracted a lot 
of interest since Fermi degeneracy was achieved by De- 
Marco and Jin pj. These systems are in the metastable 
gas phase, as three-body recombinations are rare. Most 
interestingly, the effective two-body interaction itself can 
be controlled via external magnetic fields. This makes 
it possible to study the system as it evolves from a 
dilute Fermi gas with weak attractive interactions to 
a bosonic gas of diatomic molecules. This transition 
from a superfluid BCS state to Bose Einstein conden- 
sation_(B_EC} has_been the su b jec t of many experimental 
flCA Jill jia and theoretical works 
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At the midpoint of this transition, the two-body sys- 
tem has a zero-energy bound state, and the scattering 
length diverges. If other parameters as the effective range 
of the interaction can be neglected, the interparticle spac- 
ing becomes the only relevant length scale. This defines 
the unitary limit. In this limit, the energy density is pro- 
portional of that of a free Fermi gas, the proportionality 
constant denoted by £. Close to the unitary limit, cor- 
rections are due to a finite, large scattering length a and 
a small effective range r$ of the potential. Within the 
local density approximation (LDA), the energy density 
is given as 



£[p] = £fg U 
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is the energy density of the free Fermi gas. The universal 
constant £ has been computed by several authors. Monte 
Carlo calculations by Carlson et al. [2l|. Astrakharchik 
et al. , and by Bulgac et al. [2j| agree well with each 
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other and yield £ « 0.44 ± 0.01, £ « 0.42 ± 0.01, and 
£ w 0.42, respectively. A calculation by Steele J24[ based 
on effective field theory yields £ = 4/9, while an appli- 
cation of density functional theory (DFT) [2^, |2(| yields 
£ as 0.42 [23. Other calculations deviate considerably 
from these results. Heiselberg [ljj obtained £ = 0.326, 
while Baker found £ = 0.326 and £ = 0.568 from 
different Pade approximations to Fermi gas expansions. 
Engelbrecht et al. j2£j obtained £ = 0.59 in a calculation 
based on BCS theory, while a very recent Monte Carlo 
simulation by Lee [3(j yields £ « 0.25. The experimen- 
tal values are £ « 0.74 ± 0.07 @, £ = 0.51 ± 0.04 [H, 
£ « 0.7 4], £ = 0.27j££g @. The constant c x in Eq. (JTJ 
has also been determined. The Monte Carlo results by 
Chang et al. |31l and by Astrakharchik et al. [22| yield 
ci fts —0.28 [23 an d are very close to Steele's analytical 
result [2i|. We are not aware of any estimate for the 
constant c 2 in Eq. Q that concerns the correction due 
to a small effective range. It is the purpose of this work 
to fill this gap. This is particularly interesting as exper- 
iments also have control over the effective range. Note 
that the regime of a large effective range has recently 
been discussed by Schwenk and Pethick |20j . 

In this work, we determine the coefficients c\, and c 2 
via density functional theory. Recall that the density 
functional is supposed to be universal, i.e. it can be used 
to solve the A^-fermion system for any particle number N, 
and for any external potential. Exploiting the universal- 
ity of the density functional, the parameters ci and c 2 can 
be obtained from a fit to an analytically known solution, 
i.e. the harmonically trapped two-fermion syste m [ 33ll . 
This simple approach has recently been applied |27| to 
determine the universal constant £, and will be followed 
and extended below. 

Let us briefly turn to the harmonically trapped two- 
fermion system. The wave function u(r) in the relative 
coordinate r = r\ — r 2 of the spin-singlet state is given 
in terms of the parabolic cylinder function U(—e,r/X) 
[lilli^lllll. Here, efku is the relative energy, and A = 
y/h/ (mid) denotes the oscillator length. We are dealing 
with a short-ranged two-body interaction and quantize 
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the energy through the boundary condition at the origin 



d r u(r) 



u(r) 



k cot S, 
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r=0 



where h 2 k 2 /m = eftw, and <5 denotes the s-wave phase 
shift. The evaluation of Eq. @ for the parabolic cylinder 
function yields 



V2 



r(3/4-e/2) _ A r e 
r(l/4-e/2) ~ a ~ 2A ' 
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Here, we have employed the effective range expansion of 
the phase shift. Note that Eq. Q is valid for arbitrary 
values of the scattering length a and the effective range 

As an introductory example, we consider the case of 
a dilute Fermi gas with a small value of the (positive) 
scattering length a <C A and zero range. We expand 
Eq. J3J around the energy of the noninteracting system 
as s = 3/2 + Ae. The energy correction fulfills Ae <C 1, 
and we find 



Ae 



2 a 
7T A' 
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The form of this result suggests that the energy density 
of the weakly interacting system is that of the noninter- 
acting system plus the term 



A£[p] =c[ap 1 / 3 ) -p 5 / 3 
m 



(0) 



which is due to the scattering length. We want to deter- 
mine the coefficient c in Eq. 10. Recall that Kohn-Sham 
DFT is variational, and that we are dealing with a small 
perturbation a <C p 1 / 3 . Thus, we can insert the density of 
the noninteracting system p(r) = 27r~ 3 / 2 A~ 3 e~ r / A into 
Eq. @ and integrate over all space. Equating the result 
with the energy correction given by Eq. © yields c = tt, 
which is in full agreement with many-body perturbation 
theory [H l37l l3l| . This result is not really surprising. 
The interaction is a contact interaction, and the energy 
correction given by Eq. © is the Hartree-Fock approxi- 
mation of this interaction. Nevertheless, it is encouraging 
that the simple DFT approach via the two-body system 
yields a result in agreement with many-body theory. 

Let us turn to the vicinity of the unitary regime. Con- 
sider the case of a large scattering length a » A and 
zero range. We expand Eq. 10} around the energy cor- 
responding to the unitary regime as e = 1/2 + Ae, and 
find 



Ae 



2 A 
tt a 



(7) 



This expression suggests that the correction to the energy 
density ££fg is of the form 



A£ 1 [p} = 



ci 
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FIG. 1: Energy per particle (in units of the free Fermi gas) as a 
function of (fea) -1 in the vicinity of the unitary regime. Solid 
line: slope estimated in this work; data points: Monte Carlo 
results from Ref. [2lJ (dots) and from Ref. (squares), re- 
spectively. 



We insert the exact density at the unitary regime, 
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into the correction given by Eq. © and integrate. Equat- 
ing the result with the exact result Q) yields c\ — —0.244. 
Monte Carlo calculations predict c\ « —0.28. Our result 
deviates only 13% from the results of the Monte Carlo 
calculations (see Fig.^l. The deviation is due to the fact 
that the simple functional in Eq. (JSJ is the LDA of the 
(unknown) exact density functional. Given the simplicity 
of our approach, the estimate is remarkably accurate. 

Let us consider the corrections due to a non-zero effec- 
tive range r$ <C A. Again, we expand Eq. around the 
energy of the unitary regime as e = 1/2 + Ae, and find 



Ae 



J_rp_ 
8^ A ' 
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The form of this energy correction implies that the term 

A£ 2 [p] =c 2 r p 1/3 £ F G[p] (11) 

has to be added to the energy density ££fg- For a de- 
termination of the coefficient C2, we insert the density 
given by Eq. @ in Eq. i|ll[l and integrate. Comparison 
of the result with the exact result (|10|) yields c 2 = 0.142. 
This is one of the main results of this work. We esti- 
mate that the systematic error of this coefficient is about 
5%-15%, as this is the deviation by which the DFT esti- 
mates for £ 1271 and c\ deviates from the Monte Carlo 
predictions [2l], |2^. The estimate for c 2 enables us to 
discuss a small systematic correction of the universal 
constant obtained from Monte Carlo calculations. Re- 
call that the Monte Carlo calculations [2l| and [2^] are 
based on potentials with a small effective range of about 
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rop 1 / 3 ~ 0.05, and rop 1 ^ 3 = 0.01, respectively. This sug- 
gests that their predictions for the universal constant £ 
involve a small positive error of about C2Tqp 1 ^ 3 ~ 0.007 
and C2rop 1 / 3 w 0.001, respectively, which is within the 
statistical error of these simulations. 

We also tried to improve the accuracy of our estimates 
for ci and ci by going beyond the LDA. The main idea 
consists of adding gradient terms to the energy func- 
tional, and to use Kohn-Sham DFT. The systematic in- 
clusion of the nonlocal kinetic energy density in the en- 
ergy functional can lead to improvements in the density 
and energy spectrum j^, 0] • Here, we follow a phe- 
nomenological approach. We replace the functional in 
Eq. 0) by the functional 



ci 



£[p] = + —^£a[p] + c 2 (r p i/3 )£ ro [p]- (12) 



Here 
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and similar expressions with parameters f a and f ro are 
employed for the terms involving the scattering length 
and the effective range, respectively. Note that the func- 
tional m is the Thomas-Fermi approximation of the 
functional l|12l) . and that both functionals are identical 
for = f a = f ro = 0. Note also that the density- 
dependent term in Eq. (|13fl is the Thomas-Fermi approx- 
imation of the corresponding gradient term. The pair 
of parameters (£, /f) was determined in Ref. [53, and 
the universal constant £ varies only very little when 
is varied. This is very different for the parameter pairs 
(ci,/o) and (c2,/r ), as the energy obtained from inte- 
gration of the gradient term in Eq. H13|) differs by a factor 
of 2.1 and 0.7, respectively from the energy of the corre- 
sponding density-dependent term. This finding indicates 
that the functionals £ a [p] and £ ro [p] exhibit considerable 
finite-size corrections (as the gradient terms differ from 
their respective Thomas-Fermi limits for the two-body 
system). For this reason, we do not use phenomenological 
gradient corrections for a more accurate determination of 
the constants ci and ci. 

Let us also investigate the deep bound-state limit (e — > 
— oo ) of the two-body system corresponding to a posi- 
tive scattering length a <C A and zero range. Taking this 
limit in Eq. gj), and noting that T(x + 1/2) /T(x) — > y/x 
for x — > oo, we find that the binding energy is ehuj = 
—h 2 /(ma 2 ). Thus, one can trivially write down the den- 
sity functional for the system in this limit as 



£b[p] = - 



2 ma 2 



2m a 2 



(14) 



Interestingly, this 



and the energy per particle is 

value coincides exactly with the 1/a 2 correction that Bul- 
gac and Bertsch (3^ obtained from a fit to Monte Carlo 
results close to the unitary regime, and it is about 20% 
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FIG. 2: Energy per particle for symmetric neutron mat- 
ter as a function of the density. Full line: Friedman- 
Pandharipande |4l|: dashed line: result from DFT; dotted 
curve: Fermi gas in unitary regime. The inset shows the 
DFT result and includes finite range corrections. 



larger than the analytical result that can be inferred from 
Steele's work p^ . 

Finally, we apply Eq. JQ) to neutron matter, for which 
a = —18.3 fm and ro = 2.7 fm. We drop the rop 1 ^ 3 term 
in Eq. l|T]l. as this correction is only small for very small 
densities. In Fig. [3 we compare our results to the equa- 
tion of state (EOS) by Friedman-Pandharipande [4ll |. 
Note that that EOS is based on a realistic Hamiltonian, 
which includes higher partial waves and three-body in- 
teractions. Recall that our approach is limited to s- waves 
and two-body interaction. The inset of Fig. shows the 
comparison for very small densities; here, the correction 
due to the effective range is included, and the restriction 
to s-waves is justified. We note that the inclusion of the 
effective range correction for values of r^p 1 ^ 3 less than 
0.6 improves the DFT result. 

To summarize, we have considered interacting dilute 
Fermi systems near the unitary regime and computed the 
corrections to its energy density due to a large scattering 
length and a finite effective range of the two-body inter- 
action. Our calculations are based on the universality of 
the density functional, and we determine its local density 
approximation through comparison with exact results for 
the harmonically trapped two-fermion system. The cor- 
rection due to the large scattering length agrees well with 
results from Monte Carlo calculations and effective field 
theory, while the correction due to the finite range im- 
plies a small systematic correction of order 0.01 to the 
universal constant extracted from Monte Carlo results. 
The phenomenological inclusion of gradient terms is dif- 
ficult due to finite-size corrections. We also applied our 
results to neutron matter. 
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